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APPLICATION OF KALMAN FILTERING TO ERROR CORRECTION 


OF INERTIAL NAVIGATORS 
By Heinz Erzberger 
Ames Research Center 


SUMMARY 


The design and performance characteristics of an optimum error damping 
system for an inertial navigator are investigated. The chief component of 
this system is a Ka lma n -Bucy filter which gives best estimates of the inertial 
navigator* s errors from noise-contaminated auxiliary velocity or position 
measurements. The errors estimated by this system include random and constant 
gyro drift, azimuth, leveling, accelerometer, as well as velocity and position 
errors. The system also estimates time-correlated errors in the auxiliary 
navigation measurements which are used to correct the inertial system. 

The estimates of error obtained after each measurement from the Kalman - 
Bucy filter are treated, for the purpose of error reduction of the inertial 
system, as if they represented exact error values. This procedure of using 
the estimates as if they represented perfect error measurements is also 
opt imal . 

By means of digital computer simulation, the performance of the optimum 
error damping system was compared, whenever possible, with that of conven- 
tional methods. Thus, when an auxiliary velocity is used (e.g., doppler 
radar), the optimum estimator reduces the inertial position and azimuth error 
more than 50 percent below the values obtained with the best conventional 
methods, such as velocity damping or gyrocompassing . The comparative gains 
made by use of the optimum system depend strongly on the accuracy of both the 
auxiliary measurement and the initial alinement of the inertial system. But 
it was generally found that the lower the accuracy of either auxiliary mea- 
surement or initial navigator alinement, the better is the performance of the 
optimum system in comparison to conventional methods. 

The performance of an inertial system, optimally updated with auxiliary 
position measurements, was also studied. The use of auxiliary position mea- 
surements, such as LORAN or TACAN, has not previously been exploited in a 
systematic manner; thus, no performance comparison with conventional methods 
was possible. It was found that because of the high position accuracy obtain- 
able with these navigation aids, all errors in the inertial system, including 
azimuth, were strongly reduced. In view of the relatively low cost and gen- 
eral availability of these aids, this result is felt to be especially signifi- 
cant for the projected use of inexpensive inertial systems in commercial jet 
aircraft . 


INTRODUCTION 


It is well known that an inertial navigator without periodic updating 
will eventually make unacceptably large errors. Gyro drift, mechanization, 
and platform torquing errors invariably degrade the accuracy of the system 
until, after a period of time, it ceases to be a useful navigation aid. Only 
at great cost - by careful construction of critical components such as gyros 
and accelerometers - can its period of usefulness be extended to more than a 
few hours . 

The high cost of accuracy in inertial navigators has led to the use of 
auxiliary navigation data for controlling the rate of error propagation. 

Such a system, although no longer purely inertial, often is more economical 
to construct because of reduced sensitivity to critical components. In addi- 
tion, it is more versatile than a purely inertial system since it permits at 
least partial alinement of the navigator in flight. 

The most common error-damping scheme in operational aircraft inertial 
systems is the velocity damped mode, which uses an independent velocity mea- 
surement, such as provided by a doppler radar. Other schemes also dependent 
upon some form of velocity measurement are the gyrocompassing and automatic 
leveling modes . The analysis of such damped inertial systems is treated 
elegantly by means of classical feedback control theory in reference 1. 

Although the classical theory explains the operation and aids in the 
design of the system, it is seriously deficient in some aspects. For 
instance, it does not guarantee that all auxiliary navigation measurements 
are optimally processed for error damping purposes, nor does it possess an 
efficient facility for handling measurements contaminated by noise. These 
and other difficulties with the classical method suggest taking a fresh look 
at the error damping problem in the context of modern estimation and control 
theory . 

Here the recent work of Bona and Smay on the optimum reset of ship's 
inertial navigation systems should be mentioned (ref. 2). The present paper 
differs from theirs in scope and content in that it is slanted more toward 
aircraft inertial systems . 

Essentially, the modern approach separates the error damping problem 
into two distinct operations of which the first consists of optimally esti- 
mating the errors in the inertial system from imperfect measurements, and the 
second of using the estimated errors to correct the inertial system. 

This approach has advantages both in theory and in practice. For 
instance, random processes in the error model of the inertial system, such as 
random gyro drift, for example, are handled with ease, as are noise- 
contaminated navigation measurements. Then, too, there is no limit to the 
number of sources of auxiliary noisy navigation measurements that can profit- 
ably be used for error reduction. Furthermore, the theory provides a simple 
algorithm for computing the best estimates of all the error states included 
in the error model. For instance, the algorithm gives a best estimate of 
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gyro drift rates from position measurements, such as LORAN or TACAN, even 
though this information may be available only at intermittent time instants. 
Finally, all sources of auxiliary navigation measurements are processed opti- 
mally, in the sense that the mean-square errors in estimating the inertial 
system errors are minimized. 

This paper serves a twofold purpose. First, the design of the optimum 
estimator and controller is described for a typical inertial system. Then 
the performance of optimum and classical error damping is compared by means of 
digital computer simulation. Since the two methods are compared for low and 
high accuracy inertial systems as well as for different initial alinement con- 
ditions, a rather complete picture of the characteristics of the optimum 
system is obtained. 


SYMBOLS 


Note: 

A 

C 

F 

£>m 

K± 

k 

Pi 

P i 

Qi 

Q 

R 

Ri 

u 

u e (t) 


Underlined quantities represent vectors . 
specific force vector 
measurement matrix 
system matrix 

force per unit mass due to gravity 
Kalman gain matrix 

torquer scale factor error (k x , ky, k z ) 

covariance of error states at time i + 1 given the measurement at 
time i 

covariance of error states at time i given the measurement at time 
i; also, covariance of estimation error 

covariance of noise for a discrete system 

covariance of white gaussian system noise 

variance of u £ (t) 

radius vector from center of the earth to location of navigator 
covariance of measurement noise at time i 
system noise vector 

gaussian noise generating random drift rate 
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Xd 


v 

X 


X 


11 


z 


a 

P 


V 


V 

av ; 1 

5R 

5Vr 

SVp 

50 


6 

G f 


£c 


6 

— r 


true velocity vector 

difference between inertial and reference velocities 
reference velocity vector 

random velocity reference error (v x , Vy, v z ) 
error state vector 

optimum error estimate at time i given the measurement at 
time i 

measurement vector 

inverse of correlation time of accelerometer bias error 
inverse of correlation time constant of gyro 

gain constants occurring in classical velocity damped schemes 

total accelerometer error (v x > Vyj> V z ) 

accelerometer bias error 

inertial navigator velocity error 

inertial navigator position error 

constant velocity reference error 

total reference velocity error 

vector angle relating computer set of axes to true set 
effective gyro drift rate vector Gy, g z ) 

gyro drift rate vector Gy, g z ) 

effective constant drift rate vector (including gyro torque 
errors) 

constant drift rate of a gyro 
polar component of drift rate vector 

effective random drift rate vector (including gyro torque errors) 

random drift rate of a gyro 

latitude 
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p 

2 

a 

1 

$(t,T) 

i 

a 

U) 


2c 


OJ 

-p 


Wr 


( ) 


angular rotation rate of true set of axes with respect to an earth- 
fixed. set 

steady-state variance of drift rate 

vector angle relating platform set of axes to true set 
transition matrix for error dynamics 

vector angle relating platform set of axes to computer set 

angular rotation rate of earth with respect to inertial space 

angular rotation rate of true set of axes with respect to inertial 
space 

angular rotation rate of computer coordinates with respect to 
inertial space 

angular rotation rate of platform with respect to inertial space 

Schuler angular frequency 

expected value of quantity in parentheses 


STATISTICAL MODELS FOR GYRO DRIFT RATE ADD ACCELEROMETER ERROR 


The error propagation in inertial navigation systems has been studied 
and documented by many authors (refs. 1 and 3). Several different versions 
of the error equation exist, but for the purposes of this note the particular 
form obtained by J. C. Pinson in reference 1 is most suitable. The notation 
adopted here also conforms as much as possible with that found in reference 1. 
A brief account of this theory and the notation can be found in appendix A. 
Although the error equations developed in appendix A suffice when gyro drift 
rates and accelerometer errors are deterministic time f met ions, they need to 
be augmented with additional equations of the error processes when this is 
not true . 

Consider first gyro drift rate e* , which is one of the most difficult 
error sources to control in inertial systems. In the relatively low acceler- 
ation environment assumed here the drift rate is generally composed of a con- 
stant or bias component and a slowly changing random component. To minimize 
navigation errors, it is usually necessary, prior to entering the normal 
navigation mode, to adjust the compensating torques to the gyros so that the 
total drift rates, constant plus random, are as close to zero as possible. 
Nonetheless, even with careful adjustment, drift rates of both types will 
remain. Experience indicates that random drift usually predominates and has 
the most damaging effect on the navigator accuracy. The random nature of gyro 
drift has long been recognized, but exhaustive studies of the statistical 
properties of random gyro drift are difficult to find. However, among the 
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studies available, the consensus seems to be that random drift rate of a 
single degree of freedom gyro is exponentially time-correlated (ref. 4). 
that case its autocorrelation function has the form 


G(t) 


2 - T 
= cr r e I 1 


P 


In 


( 1 ) 


The correlation time constant l/j3 lies in the neighborhood of 3 to 10 hours 
for most gyros. This is the basic statistical model for random gyro drift 
assumed here. Furthermore, if the distribution of the drift rate at any 
fixed time instant is gaussian, then the Kalman-Bucy theory applies rigorously 
(ref. 5) • A recent statistical analysis of drift data from 50 gyros would 
tend to support such an assumption (ref. 6) . The gaussian assumption is quite 
convenient, for then it is a simple matter to synthesize a linear dynamical 
system excited by white gaussian noise so that its output has the desired 
statistical properties (ref. 5)* For the random process considered here the 
proper choice of system is 


de£ 

IF 


= -Pe£ 


+ u, 


r(t) 


( 2 ) 


where denotes random drift rate. The variance, a^, of the steady-state 

drift rate can be shown to be related to the variance, q G , of the white 
gaussian noise u^t) by the equation 


d e = 2 P aJ 


( 3 ) 


Since there are three gyros on the stabilized platform, similar dynamical 
systems are assumed to generate the random drifts for each of the three 
channels . 


The constant drift rate, also known to be present, has not yet been 
discussed. Although it is not a random process, its value is nevertheless 
unknown a priori and, therefore, can only be described probabilistically. The 
distribution of the random variable will again be assumed gaussian with 
variance a 2 . The dynamic system generating this constant random variable is 
s imply 


K 

dt 


= 0 




The drift rate model for one gyro is represented in sketch (a) . 



Sketch (a) Dynamic model for gyro drift (one single-axis gyro) . 
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Depending on the characteristics of a particular inertial system, it may he 
desirable to simplify the proposed drift model. This can easily be done by 
setting either or e£ to zero. 

Finally, it is necessary to consider the contribution of torquer scale 
factor error (k x , k y , k z ) to the total drift. As shown in equation (A3), the 
scale factor error is proportional to the torquing rate w. This error is 
apparently more difficult to handle than the drift rate e f since u> is a 
time-varying quantity for a moving navigator. Nevertheless, in real time 
error estimation no difficulty is encountered. One constructs a statistical 
and dynamic model for each k, just as for the pure drift rate e* and multi- 
plies k by its respective component of w, which is available in the iner- 
tial system 1 s computer. Usually, a simple dynamic model, such as for constant 
gyro drift, should be sufficient. If w is not time varying, as in the simu- 
lation discussed later, the two sources of drift need not be considered as 
separate entities, and thus can be represented by a single model, such as 
shown in figure 1. Henceforth, gyro drift shall mean the total effective gyro 
drift e, appropriately modified to include the scale factor error. 

The above discussion on the statistical representation of total gyro 
drift applies also to total accelerometer error v- This error is com- 
posed of an offset error and a scale factor error, the latter being propor- 
tional to the specific force vector A as shown in equation (A10) . Generally 
speaking, accelerometer error is a less serious problem than gyro drift 
(ref. l) . For example, taken by itself it contributes only a bounded position 
and velocity error; f urthermore , it generally changes less with time than gyro 
drift. Thus a simple statistical model shown in sketch (b) and similar to the 



Sketch (b) Dynamic model for error of one accelerometer. 


one for random gyro drift is also chosen to represent the total accelerometer 
error v- The correlation time l/a depends, of course, on the particular 
accelerometer, but should lie somewhere between 1 and 20 hours . Strictly 
speaking, this model generates only the offset error v 1 and not the total 
error v (see eq. (A10)), except when the specific force vector A is 0. 
However, under cruise conditions and a two axes locally level platform mech- 
anization, the average value of A along the accelerometer sensitive axes 
will be small so one is justified in neglecting the scale factor error. 

One of the advantages of the Kalman filter over conventional error damp- 
ing schemes is its ability to account properly for measurement noise. In fact. 
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if the measurement noise is time correlated, then one can construct a 
dynamical model of the noise just as was done for gyro drift and estimate the 
time-correlated part of the noise along with the inertial system errors. A 
case in point occurs when doppler velocity is used for inertial system error 
damping. Doppler noise at a fixed reference velocity is composed of a ran- 
domly varying component, v, with correlation time constants in the order of 
seconds and a constant bias component SV^ . Since both noise components are 
time correlated, both could be estimated. However, if the sampling interval 
of the Kalman-Bucy filter is much longer than the correlation time of the 
noise, it is more efficient to prefilter the doppler signal with conventional 
smoothing filters. For sampling intervals of a minute or longer, it will 
still be approximately true that the random noise component of a smoothed 
doppler signal is uncorrelated. The bias component is easily included in the 
filter design if the inertial system error model is augmented with the simple 
dynamic model, developed earlier for a constant random variable, to yield the 
three equations. 


d&Vrx _ n dSV^y _ d5V r z 

“at ' “at ' “at 


(5) 


where 5V r = column (&V rx , 5V r y, 5V rz ) is related to reference velocity V r 
and true velocity V as follows: 

Xr = V + SVr + v (6) 


In effect, the Kalman-Bucy filter will now estimate the constant reference 
velocity error along with the inertial system errors. 

The next step is to write the entire set of error equations composed of 
equations ( 2 ), (4), (5), (A5), and (A9) as one first-order vector differential 
equation: 

-= = Fx + u ( 7 ) 

dt 


The ordering of the error variables into a state vector x is, of course, 
arbitrary but a logical choice consists of the following arrangement: 


x = column 

t G rx> 


£ rz> 

€ CX^ 

e cy> e cz> V *z> 

V X , VyJ 

AV X , 

AVy, 

SR*, 

SRy, 

6 V rx , SV ry , &V rz ] 


( 8 ) 


Equation ( 8 ) implies the definition AV X = &R X , AVy - 5R y and the assumption 
that the vertical channel is not mechanized, that is, 5R Z = 0. The matrix F 
is, therefore, of dimension l 8 xi 8 , and its entries are easily determined from 
the various error equations. Finally, the noise vector u has the form 

u = column [u ex , u ey , u ez , 0 ,0,0, 0,0,0, u 7X , Uyy, 0,0, 0,0, 0,0,0] ( 9 ) 
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MEASUREMENT EQUATIONS 


l| 
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The purpose of this section is to develop the relationships between the 
measured error and the error state variables. Such relationships are essen- 
tial components of the Kalman -Bucy filter theory. 

Consider first an external velocity source , such as doppler radar , used 
as a reference . The doppler radar gives the components of the aircraft veloc- 
ity vector with respect to the earth in a coordinate frame fixed in the air- 
craft . The components of this velocity vector are transformed to platform 
coordinates with the help of the platform resolvers which specify the relative 
orientation of the two coordinate systems. Since platform and computer frames 
are nominally identical, inertial and doppler velocity can now be compared. 

In general, the resulting difference velocity contains elements of 
inertial as well as doppler velocity error. The expression for the error of 
inertial system velocity with respect to an earth-fixed frame can be obtained 
directly from the equation for V derived in reference 1: 

V = rirRl +P XR 

— Ldt Jc ~ c ~~ 

where is the rotation rate of any given set of axes with respect to an 

earth-fixed set, and the time derivative of R is also taken with respect to 
the given axes. The error in V, denoted by AV* , is then obtained from 

equation (10) by substituting R + SR in place of R and identifying terms: 

V + AV‘ = [*— 6r] + p X 6R + [— r] + p X R 

— Ldt — J ~c — Ldt - I -c - 

C c 

= ZW+£X5R+R + pXR (11) 

The right side of equation (ll), resolved along some known set of axes, 
usually the computer set, generates the actual velocity available in the com- 
puter. This inertial velocity is to be compared with the reference velocity 
Xr = X “ ^Xr viiere SV j. ~ &Vr + X denotes the total error in the reference 
velocity. But the existence of an error angle jr between computer and plat- 
form axes means that the computer will actually utilize 

V - 5V^ - t X (V + 8V*) « V - 6V£ - jf X V 

in the velocity comparison. Therefore, the difference between inertial 

and reference velocities is 

V d = ^ + pX5B+tXV + 5V^ 1 (12) 

In the preceding section the reference velocity error was separated into bias 
and white noise components and the bias component was included in the inertial 
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system error model. If this is also done in equation (12), one obtains the 
final form of the measurement equation 

V d =AV;+£X 5R+tXV + 5V r + v (13) 

■where v denotes the time -uncorrelated (white noise) component of reference 
velocity error. 1 Since equation (13) expresses the measured quantity as 

a sum of uncorrelated noise and a linear combination of the error state vari- 
ables, it can be put in the form required for the Kalman-Bucy theory (ref. 4); 
that is to say, a 3x18 matrix C can be found which allows equation (13) to 
be written as 


z = Cx + v (l4) 

Here x represents the error state vector and z the measured quantity, in 
this case . 

Another important reference source consists of auxiliary position infor- 
mation such as might be obtained from LORA N, TACAN, or even from visual obser- 
vation of known landmarks. In the past these measurements, although readily 
available in many aircraft, have not been fully utilized for error damping of 
inertial navigators. Yet accurate position information can help considerably 
in decreasing heading error, as the computer study to be discussed later 
demonstrates . 

The measurement equation for position reference is simply 

% = 8R + v p (15) 

The vector Vp represents the reference position error assumed to be white 
gaussian noise. Obviously, equation (15) can also be written in the form of 
equation (l4) . 


DESIGN OF THE FILTER 


Once the linear dynamic system and the matrix equation relating the state 
variables of the system with the measured quantities have been determined, the 
Kalman-Bucy filter can be constructed, provided, of course, that certain sta- 
tistics, such as the covariances of the system noise and measurement noise and 
the initial covariance of the state variables, are also known. However, 
before the filter equations can actually be written, two questions must still 
be resolved. The first is how the existing filter theory can be made appli- 
cable to continuous systems sampled at discrete time instants; and the second 
is how to use the estimated error states to the best advantage in reducing the 
inertial system errors. 


■*The components of v may be correlated with each other, depending on 
the particular hardware utilized. 
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The significance of the first question above is best understood by- 
envisioning the sequence of events as they would occur in the system. Suppose 
a comparison has just been made between inertial velocity and an auxiliary 
navigation measurement, say doppler velocity. The difference velocity, at 
whatever time instant it was measured, must be accepted by the filter which 
must operate upon it and yield a best estimate of the error state. Setting 
aside for the moment the question what to do with the computed estimate, con- 
sider what may happen next at the filter input. Perhaps the aircraft is over 
water and the doppler signal begins to fade, or in case of position measure- 
ment by LORAN, it has just passed outside the range of LORAN stations. In 
other words, neither the arrival time nor the type of the next measurement is 
known beforehand. Even if auxiliary navigation measurements were available 
continuously without interruption, the filter still could accept signals at 
discrete time instants only, because of the limited computing speed of the 
digital device which mechanizes the filter equations. The solution to all 
these difficulties lies, of course, in changing the time-continuous error 
equations to a discrete system, at least for the purpose of filter design, 
while at the same time taking care to modify the statistics correctly. 


The time discrete system is obtained by expressing the solution of equa- 
tion ( 7 ) in terms of the fundamental matrix $(t i+1 ,tj_) of the system (ref. 7 ) 

(16) 


x(t 1+1 ) = 0(t 1+1 ,t i )x(t 1 ) + 


I 


$( t i+ 1 > T )u( T ) dT 


This equation is written in abbreviated form as follows: 

2i+i = ®x2£i + Hi 


(17) 


To complete the transition from the continuous to the discrete system, it is 
necessary to calculate the covariance of the random vector uq in terms of 
the known system parameters and the covariance of u(t) • The covariance Qq 
of Uj_ is calculated as follows: 


^ t -? j..? 


n u i+i p 

Q ± = J J ^(t i+1 ,T)u(r)u T (t*)0 T (t i+1 ,t I )dT dt* 


*1 


t i 


(18) 


But u(t) is a white gaussian random process; hence, its autocorrelation 
matrix Q,(T,t*) has the form 


u(T)u T (t') = Q(T,t') = Q(t ' ) & (t - t’) 


( 19 ) 


where b(r - t 1 ) is the delta function (ref. 7)* Then using equation ( 19 ) and 
a property of the delta function simplifies equation (l8) to 


Q 


1=1 $(t i+1 ,T)Q(T,T) l I> T (t i+1 ,T)dT 

*i 


( 20 ) 


which gives the desired relation. 
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With these modifications it is clear that the Ralman-Bucy filter theory 
for discrete systems is directly applicable here. Since the construction of 
the filter and the form of the recursion relations for computing the filter 
parameters are readily available (e.g., refs. 8 and 9 ) 9 only the final results 
are included here. Let Xjl ^ denote the optimum estimate of Xi given the 
measurement z±. The next estimate x^+i| i+i is computed from the previous 
estimate x^j ^ by means of the following relation: 

2Ei+i|i+i = + K i+1 [z 1+1 - (21) 

Here C^ +1 is the measurement matrix, which is interpreted in equation (21) 
as the value C assumes at the (i + l)st time instant . The matrix 
is called the gain matrix of the filter and is computed recursively as follows: 

P ± = + Q ± (22) 

P± +1 (1 ^i+i^i+i^l (23) 

Ki+i " P i C i +1 <°i«. P i C l+l + R 1+1 ) _1 (24) 

where 

P^ = (x.j_ - & i)i )(x i - 5c_j_ I j_ ) T (covariance of the estimation 

‘ 1 error at time i) 


T 

Ri +1 = Xi-fiZi+i (covariance of zero mean 

measurement noise at time i + l) 

Pi is defined by the right side of equation (22) . To start the itera- 
process at i = 0, one must have available the a priori statistics P Q 

x o|o* 

Now consider the answer to the second question posed at the beginning of 
section, namely, how to use the error estimates to the best advantage in 
reducing the inertial system errors. Since the optimum estimates represent, 
in a certain sense, the best information available on the error states up to 
that time, one might argue that the optimum decision after an estimate has 
been computed consists of resetting every inertial system variable by the 
amount of the error estimate. Specifically, this procedure would entail sub- 
tracting 6R X from the inertial x coordinate of position immediately after 
5R X has been computed and treating the other variables in a similar manner. 
That this is indeed the optimum strategy to pursue follows directly from 
published results on the combined estimation and control problem (ref. 9 ) • 

In other words, when formulating the error control of an Inertial system as a 
stochastic control problem with the cost function quadratic in the error 
states, one finds that the solution consists of an independently designed 
Kalman-Bucy filter for estimating the error states and a controller that oper- 
ates upon the estimates as if they represented perfect measurements . 


and 

tion 

and 


this 
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Resetting the inertial system variables closes the loop around the 
Kalman-Bucy filter, all the estimated error states being reduced to zero when- 
ever data are processed by the computer. This procedure simplifies the filter 
structure because at each step the first and last terms in equation (2l) are 
zero: 


— i+i| i+x K i+x^i+i 

Equation (25) holds even if the resetting does not take place instantaneously, 
so long as it is completed before the next measurement is accepted by the 
filter . 

There is, however, one complication associated with the proposed reset 
procedure that should be mentioned. It will be remembered that random gyro 
drift rate e r and accelerometer bias error v were modeled by dynamic sys- 
tems of the type shown in sketches (a) and (b) . These systems do not actually 
exist as physical entities; and therefore, it is not possible to reset them 
physically (i.e., by the conventional method of applying an appropriate con- 
trol signal at the input of the integrator) . The resetting can nonetheless be 
carried out indirectly by adding a signal to the output variable , which does 
have physical existence in the inertial navigator. The necessary reset signal 
is a slowly time-varying quantity dependent upon the dynamical properties of 
the system. Moreover, it is easy to see that the required signal is generated 
by a system identical to the one to be reset. For the system shown in 
sketch (b), for example, the signal which must be added to v is ve -<x *k if 
an observer at the output is to measure the same response as at the output of 
a system whose state had been increased instantaneously by an amount 0 at 
time t = 0. This concept is illustrated in sketch (c) below. Obviously the 
same considerations apply also to the gyro drift model. 




The design of the filter and its relationship with peripheral systems is 
depicted in figures 1 and 2. Figure 1 shows diagrammat ically the sequential 
computation of the optimum estimate, while figure 2 shows the complete closed 
loop system, including the simplified calculation of the estimate resulting 
from the procedure of setting the estimated error state to zero after each 
estimation. Not shown in the figures are the several logical operations 
necessary to integrate the filter into the inertial system as a whole. These 
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operations consist of sensing the type of measurement that has occurred, 
whether LORAN, TACAN, doppler velocity, or other, programming the computer to 
calculate the correct gain matrix K , and storing the- measurement until the 
gain matrix has been calculated. 

Not only are the parameters of equation (17) and of the measurement 
matrix C time varying, but the latter may even change dimension randomly at 
the measurement time. The calculation of the gain matrix K± must, there- 
fore, be done in real time. For the particular kind of problems considered 
here, it was found that the basic iteration interval should be chosen in the 
order of a few minutes or less so that the fundamental matrix can be 

approximated by a first degree Taylor series expansion of e F ( ti+i-tj T Also, 
over such short time intervals it is sufficient to take F constant and 
update it only at each iteration. 


EXPERIMENTAL RESULTS 


In this section inertial system error estimation by means of a Kalman - 
Bucy filter is compared with conventional error damping methods by means of 
machine computations. The results should help a designer in judging how 
greatly performance is improved with a Kalman-Bucy filter over conventional 
methods. Particular care was taken to compare the modern and classical 
methods under the same conditions so as to assure a fair comparison of the 
performance characteristics. The machine computations, summarized here in the 
form of graphs, are essentially time histories of inertial system rms error 
states calculated for a variety of operating conditions. A conventional 
locally level, latitude -longitude mechanization was chosen as the model for 
the inertial system. 

Two types of reference sources were considered, doppler velocity and 
position information. With doppler velocity as a reference, the error state 
rms values, given by the square root of the diagonal elements of P^ , were 
calculated for the Kalman-Bucy filter configuration, for velocity damping, 
automatic leveling and gyro compassing taken together, and for velocity damp- 
ing alone. Since there are no conventional error damping methods dependent 
upon position information except for elementary position reset, only the opti- 
mum estimator performance was calculated when position information was used. 

An iteration interval (tj_+i - t^) = 2*5 minutes was used for all computations. 
Also, it was assumed in this study that the vertical channel is not mechanized, 
so that the position error equation (A 9 ) resolved along the computer set of 
axes reduces to two scalar equations: 

^x^ ^s I — I — i— I ] &R x — 2co2;SRy 

+ [~p z + (uiy + Q y )p x ]6R y = 

BRy + [( w y + fty)p y + W 1 + | | “ | w 2 | ] &Ry + 2w z &R x 

[ Pz ^x ) Py ) SRx = 


^Z-^y ^yA-z Vx 


> (26) 


^X A Z “ ^Z A X + Vy 


Ik 



For doppler velocity reference, equation ( 13 ) resolved along the computer set 
of axes yields three scalar equations: 

Vdx = AV x ~ P z BR y + t y V z - t z V y + 5V rx + v x 


V dy ~ AV y + P z 6R x + ^z V x " Vz + 5V ry + v y } 
^dz = Px^ R y “ Py^ R x + ^x^y " ^y V x + B ^rz + v z 


(27) 


Equations (A5) and ( 2 6 ) together with the dynamic models for gyro drift, 
accelerometer bias, and constant reference velocity error constitute the sys- 
tem equation ( 7 )> whereas equation ( 27 ) determines the measurement equa- 
tion (l4) * Recursion relations (22)-(24) along with the system and 
measurement equations were employed in calculating the time histories of the 
rms estimation errors for the optimum estimator. It should be observed that 
the reset signals do not modify the covariance of the estimation error since 
they are deterministic (no error in the reset implementation) and, therefore, 
do not contribute to uncertainties within the system. 


The situation is somewhat different for conventional error damping 
schemes in that both the system equation ( 7 ) and the recursion relations (22)- 
(24) must be appropriately modified. Thus, in the velocity damped mechaniza- 
tion 7 x Xd add e d to the left side of the pure inertial mechanization 

equation (A 6 ), where the gain constant 7 , is selected to give the desired 
damping of the transient response (ref. 1) . This, in turn, modifies the 
position error equation (A 9 ) by adding 7-jVd to its left side. In the auto- 
matic leveling mechanization the gyros are precessed at rates w x - (? 2 / R ) V dy' 
Wy + ( 7 2 /R)Vax> where again 7 2 determines the transient response 

(ref. l) . In this case (~ 7 2 /R)Vay and (7 2 /R)Vdx must be added to the right 
side of the first and second of equations (A5) > respectively. Finally the 
earth-rate gyrocompass mode is obtained by precessing the z gyro at the rate 
u ) 2 + ( 7 z/R)Vdy and the x and y gyros as in the automatic leveling mode. 

The resulting modification to the \| r equation (A5) is obvious. To compute 
the covariance of the inertial system error of these configurations, it is 
only necessary to set the measurement matrix C equal to zero and then use 
equations (22)-(24) as before. 

Two calculated quantities that occur in the figures are platform azimuth 
error cp z and polar component of drift rate Gp . Both are linearly related 
to the chosen error state variables , 

5R* 

'Pz = tz + -g 5 tan 9 


Gp = Gy cos 0 + e z sin 0 0 = latitude of navigator 

The polar component of drift is important because it produces unbounded 
components (ref. l) . 
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In order to check the effect of truncation errors on the calculated 
values of the K and P 1 matrices, two FORTRAN programs were written, one in 
single precision, the other in double precision. Comparison of the matrices 
calculated by the two programs showed differences at most in the third signifi- 
cant figure after 120 iterations. Single precision should, therefore, suffice 
to mechanize the iteration equations. 

The results of the simulation now to be described are given for three 
cases: The first corresponds to a poorly alined inertial system equipped with 

high quality gyros, the second to a well-alined system of the same configura- 
tion as used in case 1, and the third case corresponds to a well-alined iner- 
tial system equipped with relatively high drift rate gyros . The performance 
of each case was calculated when operated with competing error damping schemes; 
namely, optimum filtering and classical velocity damping and gyro compassing. 
The complete list of parameters specifying each case can be found in tables I- 
III. 


The performance of case 1 is shown in figures 3 and 4 in the form of time 
histories of rms latitude and azimuth error propagation. Longitude error is 
not shown, but for this configuration did not differ significantly from lati- 
tude error. As is evident from these figures, the performance of the system 
equipped with an optimum filter surpasses that of the conventionally damped 
system in both the transient and steady-state phases of the response but the 
improvement is particularly evident in the transient phase. The strong exci- 
tation of the Schuler loop occurring under poor alinement conditions is 
responsible for the large peaked transients in both figures in the gyrocompass 
mode (ref. l) . Although this transient in the gyrocompass mode can be reduced 
by proper choice of the feedback gains y- 7 2 j 7 z; can only be reduced at 
the detriment of steady-state error in the presence of measurement noise. As 
the figures demonstrate, this basic shortcoming of the gyrocompass mode is 
completely overcome with the optimum estimator. Also shown in these figures 
are the responses of the optimally damped system when updated with position 
measurements at every iteration step (i.e., at 2*5 min intervals). Azimuth 
alinement accuracy of this configuration shows the characteristic lag that 
occurs whenever the estimated quantity is computed from derivatives of 
measured quantities. 

The performance of case 2, the same inertial system as case 1 except that 
it was assumed more accurately alined initially, is shown in figures 5 and 6. 
Clearly, the gyrocompass mode was misused in this case, because as seen in 
figure 6, the azimuth alinement accuracy of the pure inertial system was 
apparently superior to the steady-state gyrocompass accuracy. While the 
inflight usefulness of the gyrocompass mode is thus restricted, the optimum 
estimator does not exhibit this disadvantage since the relative accuracy of 
the inertial system and the measured data is automatically reflected in the 
weighting matrix K ± . It is also evident from the figures that the perfor- 
mance improvement gained with an optimum estimator over the velocity damped 
mode when the inertial system is accurately alined initially is not nearly so 
striking as in the poorly alined case. Thus, the greatest performance gains 
can be expected whenever the initial uncertainties of the error states of the 
inertial system are large. This finding confirms one’s belief that statisti- 
cal analysis increases in value as the disorder in the system increases. 
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Figures 7-9 show the performance of case 3 which corresponds to an 
inertial system moving south from north latitude 37-2° with a constant speed of 
1700 knots . The rather high gyro-drift rates of this system, given in 
table III, were chosen in order to evaluate the performance of the estimator 
when operated with a low-accuracy system. Again, the performance of the 
estimator surpasses the steady state and transient performance of the conven- 
tional velocity damped system by 50 percent or more. Of particular interest 
is the behavior of the position updated inertial system such as would be pro- 
vided by LORAN or the proposed navigation satellites. It is interesting to 
note the high azimuth accuracy attainable with a low accuracy inertial system 
when high quality position information is available. In comparison to doppler 
radar, LORAN is a relatively cheap navigation aid, is available over many 
important air routes and, when used to update the inertial system, can trans- 
form an inexpensive, low accuracy inertial system into one of good accuracy. 

Figure 9 demonstrates the ability of the estimator to estimate constant 
gyro drift rate along the polar axis from position or velocity measurements 
while in flight. The effective constant polar axis drift rate with position 
updating is reduced by as much as 25 percent from the value of the unaided 
system. A proportionally larger reduction in drift rate would be obtained 
with larger initial drift rates. Random drift rate (not shown) is also 
reduced by about the same percentage. In the past, gyro drift could be mea- 
sured during flight only with the aid of star trackers . 


CONCLUDING REMARKS 


The accuracy and versatility of an inertial system can be greatly 
improved if inertial system errors are optimally estimated with auxiliary 
noisy navigation measurements. While it is true that the estimator increases 
the complexity of the inertial system's computer, it is believed that the cost 
of greater computer complexity can be more than offset because cheaper iner- 
tial system components can be used for a specified navigator accuracy. Since 
an inertial system equipped with an optimum filter estimates all time- 
correlated errors included in the model, even gyro drift rates and accelerom- 
eter bias, from auxiliary navigation measurements, the alinement of the iner- 
tial system takes place automatically, on the ground or in the air. The 
optimally damped system is also far superior to the conventional velocity 
damped and automatic alinement modes, especially so when the measurements are 
inaccurate and the initial alinement of the inertial system is poor. 


Ames Research Center 

National Aeronautics and Space Administration 
Moffett Field, Calif., Nov. 29 j 19^6 
125 - 17 - 03 - 02-21 
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APPENDIX A 


INERTIAL SYSTEM ERROR MODEL 


This appendix provides a brief review of the inertial system error 
equations and the notation as used in this paper. A complete derivation can 
he found in reference 1. 

Consider a typical inertial system consisting of a gyro-stabilized plat- 
form, three mutually perpendicular accelerometers mounted on the platform, and 
a computer to solve the mechanization equations of the system. Three coordi- 
nate systems enter into the error analysis. 

A A A 

l X p,lyp,l Z p unit vectors along the sensitive axes of the accelerometers; 
also referred to as the platform coordinates 

A A A 

l x ,ly,l z unit vectors representing ideal alinement of the inertial sys- 

tem, the true coordinates at the navigator's location 

AAA 

l xc AycAze computer coordinates 

In addition, the following quantities are defined: 

(£ angular rotation rate of true set of axes with respect to inertial space 

i£ c angular rotation rate of computer coordinates with respect to inertial 

space 

ojp angular rotation rate of platform with respect to inertial space 

p angular rotation rate of the true set of axes with respect to an earth 

fixed set 

59 vector angle relating computer set of axes to true set 

cp vector angle relating platform set of axes to true set 

\jr vector angle relating platform set of axes to computer set 

Clearly the last three quantities are related by the equation 

= t + 69 (Al) 

The use of the term vector angle implies that only small angles are of inter- 
est; hence, the theory of infinitesimal rotation applies. As is shown in 
detail in reference 1, the error equations logically fall into two groups, 
namely, the j/ equations and the position error equations. The \| / equation 
is known to have the following elementary form: 
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e 


(A2) 



where the subscript I indicates that the time derivative of jr is to be 
taken with respect to an inertial coordinate system. The effective drift rate 
of the gyros , e , is defined as the sum of the actual drift rate, e f , and the 
drift rates contributed by gyro torquer scale factor errors, k x , ky, k z : 

£ = + + ly-^y^y + ^z-^z^z (A3) 


Here to x , Wy, to z represent the components of w in the true coordinate sys- 
tem. The value of <£, which is approximated in the computer by u) c , depends 
upon the particular mechanization as well as the speed and location of the 
navigator. An especially useful form of equation (A2) is obtained by taking 
the derivative in a coordinate system rotating with rate w with respect to 
inertial space, 

d\|r 1 , . v 

+ u) X \|r = € (A4) 

dt t “ - 


and then resolving this 


equation along the true coordinate axes: 
^x + ~ ^z^y = G x ' 

^y + ~ w x^Z ~ € y i 

+ w x^y " Vx = e z , 


(A5) 


This is the form of the 


\|r equation as used here. 


The derivation of the position error equations starts first with the 
writing of the mechanization equation for the particular inertial system that 
is being studied. This equation is simply Newton's second law written for a 
rotating coordinate system, and it has the following form (ref. l): 

r* r 1 + 2 ^c x 

Ldt 2 He 

(a6) 


+ 


X R + w X (w X R) - fl X (fl X R) 


where 

R radius vector from the center of the earth to the location of the 
navigator 

ft angular rotation rate of the earth with respect to inertial space 

A specific force vector measurable with the accelerometers 

g £m “ £ x (£ x R) > and is the force per unit mass due to gravity 
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A sufficiently accurate approximation of g for the purpose of this paper is 
given by 


where 


& = - 


R 

R “ 


g = acceleration of gravity 


- W S 3 


(AT) 


OJ 


S 



is the familiar Schuler angular frequency. Errors are introduced in mechaniz- 
ing equation (A 6) first because of the error v in the accelerometer outputs , 
and second because the accelerometers measure acceleration along the platform 
axes, whereas the computer treats the outputs as though they were measured 
along the computer axes . The combination of these two effects results in an 
actual computer input of 

A - x A + v (A8) 

instead of the exact value A. If equations (A7) and (A8) are substituted 
into equation (A 6) along with R + 6R in place of R, one obtains the posi- 
tion error equation by identifying terms: 

6R 4- 2w X SR 4- jd X SR 4- [ ((£ 4- Q) • SR ] p 4- + |n| — | w| ) SR = -J/ X A + V 

(A 9) 


Equation (A 9 ) holds for an observer not rotating with respect to the true 
coordinate system; for any other rotating coordinate system it is only neces- 
sary to replace w with the angular rotation rate, with respect to inertial 
space, of the new coordinate system. 

In a manner similar to gyro drift, one separates the accelerometer error 
into two components, an offset error, v% resulting from incomplete nulling 
under zero acceleration, and a component resulting from accelerometer scale 
factor errors k x> k z* 

V =: V* + lx k X^-X lykyAy + lz^2;A z (AID) 

Equations (A4) and (A9) specify the dynamics of error propagation in an 
inertial system, and they have the advantage over other formulations in that 
the angle and position error equations are uncoupled. A possibly important 
error source, which has not been included, is the uncertainty in the calcu- 
lated Schuler angular frequency w s . This error, introduced mainly by alti- 
tude error, is the origin of the vertical channel instability (ref. l) . 
However, here it is assumed that the vertical channel is not mechanized and 
that altitude is measured with negligible error by auxiliary equipment. 

Another point to be mentioned is that the analysis presented here gives rise 
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to the linear error equations whereas a more detailed treatment would show the 
existence of second degree terms in the error state variables. However, for 
the relatively small errors and short time intervals of interest here the 
second degree terms do not contribute signif icantly to the error propagation. 


21 



REFERENCES 


1. Pinson, John C.: Inertial Guidance for Cruise Vehicles. Ch. 4 of Guid- 

ance and Control of Aerospace Vehicles, C. T. Leondes, ed., McGraw-Hill 
Book Co • , Inc . , 1963 > PP* 113-187 • 

2. Bona, B. E.j and Smay, R. J.: Optimum Reset of Ship's Inertial Navigation 

System . IEEE Trans . Aerospace and Electronic Systems , vol . 2 , no . 4 , 
July 1966 , pp . 409-414 • 

3« Broxmeyer, Charles: Inertial Navigation Systems. McGraw-Hill Book Co., 

Inc . , 1964 . 

4. Dushman, Allan: On Gyro Drift Models and Their Evaluation. IRE Trans. 

Aerospace and Navigational Electronics, vol. ANE-9, no. 4, Dec. 1962, 
pp . 230-234 . 

5. Kalman, R. E . ; and Bucy, R. S.: New Results in Linear Filtering and 

Prediction Theory. Trans. ASME, Series D, J. Basic Engr., vol. 83, 
no . 1 , March 1961 , pp . 95-108 . 

6. Cooper, John R.: A Statistical Analysis of Gyro Drift Test Data. Master 

Thesis, MIT, Sept. 1965* 

7. Zadeh, L. A.; and Desoer, Charles A.: Linear System Theory. McGraw-Hill 

Book Co., Inc., 1963* 

8. Kalman, R. E.: A New Approach to Linear Filtering and Prediction Problems. 

Trans. ASME, Series D, J. Basic Engr., vol. 82, pt. 2, i960, pp. 35-45* 

9. Joseph, P. D.; and Tou, J. T.: On Linear Control Theory. Trans. AIEE, 

pt. II, 80, 1961, p. 18. 


22 



TABLE I.- ERROR SCHEDULE AMD SYSTEM PARAMETERS CORRESPONDING TO CASE 1 


Navigator location: 37-2° north latitude; gyro drift rate time constant 

l/P = 5 hours; accelerometer bias time constant l/a = 20 hours; observation 
time interval =2.5 minutes for optimum estimators and continuous for conven- 
tional modes; navigator stationary or moving slowly. 

The covariance matrix of the error states P' at t = 0 was assumed to 
be diagonal and its rms__ diagonal elements are listed below; covariance 
matrix of system noise Q is diagonal and its nonzero elements are computed 
from the steady-state drift rates and accelerometer bias as 2o§P and 2o^a, 
respectively. 


Error state and units 

Element of P ' 

Initial rms 

x gyro random drift rate £xr> min. of arc/hr 

P*(1,U 1 

1.55 

y gyro random drift rate £yr> mi* 1 arc/hr 

P' (2,2) 

1.55 

z gyro random drift rate e zr , min of arc/hr 

P* (3,3) 

1.55 

x gyro constant drift rate e X c> min of arc/hr 

p'(M-) 

1.55 

y gyro constant drift rate 6yc> min of arc/hr 

P' (5,5) 

1.55 

z gyro constant drift rate e zc , min of arc/hr 

P' (6,6) 

1.55 

\| r x , min of arc 

R' (7,7) 

30 

tyj min of arc 

P' (8,8) 

30 

min of arc 

P' (9,9) 

60 

x accelerometer random bias v x > knots/hr 

P' (10,10) 

20 

y accelerometer random bias Vy> knots/hr 

P* (H,ll) 

20 

x velocity (latitude) error AV X , knots 

P' (12,12) 

1 

y velocity (longitude) error AVy, knots 

P* (13A3) 

1 

x position (longitude) error 5R X , nautical miles 

P’(l4,l4) 

1 

y position (latitude) error 6Ry, nautical miles 

P' (15,15) 

1 

constant x reference velocity error, knots 

P* (16,16) 

2 

constant y reference velocity error, knots 

P» (17,17) 

2 

constant z reference velocity error , knots 

P* (18,18) 

2 


Velocity measurement noise covariance matrix 


Ri = 


0.25 

0 

0 


0 

0.25 

0 


0 

0 

0.25 


knots 


2 


Position measurement noise covariance matrix 


Ri = 


0.25 

_0 


0 

0.25_ 


nautical miles 


2 


Steady state random gyro drift rate (each channel): = 2 min/hr 

Steady state random accelerometer bias (each channel): cr-^ = 20 knots/hr 
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TABLE II.- ERROR SCHEDULE AND SYSTEM PARAMETERS CORRESPONDING TO CASE 2 
Only changes from table I are listed. 

€ zr e xc e yc e zc 'I'x 7y 'I'z 

0 T 707 0.707 O .707 0.707 0.5 0.5 1 



TABLE III.- ERROR SCHEDULE AND SYSTEM PARAMETERS CORRESPONDING TO CASE 3 


Navigator moves southward from 37*2° north latitude at 1700 knots . Only 
changes from table I are listed. 



Velocity measurement noise covariance matrix 
26 0 0~ 

Rj_ = 0 26 0 , knots 2 

_0 0 26 _ 


2k 



-I 


* -A 


Kj b>(?H 


DELAY 



INERTIAL SYSTEM 
ERROR MODEL 


Figure 1.- Schematic representation of the Kalman-Bucy filter. 


SYSTEM 

NOISE 



Figure 2.- Block diagram of optimum system. 











min of arc I DISTANCE, nm 



3.- Propagation of rms latitude (Y axis) error in a poorly alined 
inertial navigator (case 1). 
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Figure 4-- Propagation of rms azimuth error in a poorly alined inertial 

navigator (case 1). 




Figure 5 • - Propagation of rms longitude error in a well alined inertial 

navigator (case 2). 



Figure 6.- Propagation of rms azimuth error in a well alined inertial 

navigator (case 2). 
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min of arc c DISTANCE, nm 



re 7 . - Propagation of rms latitude error in a fast moving inertial 

navigator (case 3)* 



min 
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Figure 8.- Propagation of rms azimuth error in a fast moving inertial 

navigator (case 3)- 



PURE INERTIAL MODE. 

OPTIMUM ESTIMATOR WITH VELOCITY REFERENCE 

OPTIMUM ESTIMATOR WITH POSITION REFERENCE 



Figure 9 -- Constant drift rate about rms polar axis in a fast moving 

inertial navigator (case 3)* 


NASA- Langley, 1967 
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